function [lat, lon] = XYtoGPS(x, y, ref_lat, ref_lon)

    CONSTANTS_RADIUS_OF_EARTH = 6371000.;

    x_rad = x/CONSTANTS_RADIUS_OF_EARTH;
    y_rad = y/CONSTANTS_RADIUS_OF_EARTH;
    c = sqrt(x_rad * x_rad + y_rad * y_rad);

    ref_lat_rad = ref_lat*pi/180;
    ref_lon_rad = ref_lon*pi/180;

    ref_sin_lat = sin(ref_lat_rad);
    ref_cos_lat = cos(ref_lat_rad);

    if abs(c) > 0
        sin_c = sin(c);
        cos_c = cos(c);

        lat_rad = asin(cos_c * ref_sin_lat + (x_rad * sin_c * ref_cos_lat) / c);
        lon_rad = (ref_lon_rad + atan2(y_rad * sin_c, c * ref_cos_lat * cos_c - x_rad * ref_sin_lat * sin_c));

        lat = lat_rad*180/pi;
        lon = lon_rad*180/pi;

    else
        lat = ref_lat;
        lon = ref_lon;
    end